Rare events and breakdown of simple scaling in the Abelian sandpile 
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tion of SOC |11|. For ID sandpiles evidence has been 
most often definitely in favor of muftiscaling [pl]-p^. 
However, the physical origin of such multiscaling was 
never fully elucidated. On the other hand, for systems 
in rf > 1 standard FSS has never been put into serious 
doubt, so far, and, with almost no exception |T^, has 
always been assumed in both numerical and theoretical 
work. 

We consider a square lattice box of side L. At each site 
i an integer variable Zi = 1,2,... represents the number 



of 



grams. 



If 



4 the site topples, i.e. Zj 
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Due to intermittency and conservation, the Abelian sandpile in 2D obeys multifractal, rather than 
finite size scaling. In the thermodynamic limit, a vanishingly small fraction of large avalanches 
dominates the statistics and a constant gap scaling is recovered in higher moments of the toppling 
distribution. Thus, rare events shape most of the scaling pattern and preserve a meaning for effective 
exponents, which can be determined on the basis of numerical and exact results. 
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The sandpile model was introduced as a first example 
of self-organized criticality (SOC) SOC can manifest 
itself in the stationary state of slowly driven, dissipative 
systems: intermittent activity bursts occur then over all 
allowed scales, implying power law correlations in both 
space and time. 

The analytic tractability of the Abelian sandpile 
(ASM) has allowed a number of exact results, especially 
in 2D This considerably enhanced the theoreti- 

cal interest of this prototype model. On the other hand, 
in spite of some remarkable advances |§,||,||,D , a satisfac- 
tory understanding of ASM scaling is still missing. Unlike 
many stationary state properties, the exponents describ- 
ing avalanche size distributions are not known exactly, 
yet, and remain a major challenge in the whole field of 
non-equilibrium critical dynamics. 

By analyzing the avalanches of a 2D sandpile in terms 
of waves , Priezzhev et al. Q proposed a scaling picture 
which led to conjecture all exponents, including those 
describing the avalanche distributions in terms of cov- 
ered area and toppling number. Unfortunately, determi- 
nations of such exponents are notoriously difficult and 
there is no compelling numerical support of the conjec- 
ture 1^,^. A more important objection came from an 
extensive study pO[ |, which showed the illegitimacy of a 
basic assumption made in Ref. 7, concerning the relative 
sizes and scaling of successive waves. 

In this paper we reconsider scaling in the 2D ASM. 
We show that intermittency and conservation lead to a 
breakdown of finite size scaling (FSS) in the distribu- 
tion of topplings. This distribution has manifest multi- 
scaling properties. However, due to a very peculiar role 
played by a class of rare, large avalanches, a standard 
FSS description can be effectively recovered, if one fo- 
cuses on the higher moments of the distribution. Quite 
remarkably, within our picture the effective scaling expo- 
nents of the toppling number distribution are determined 
by an asymptotically vanishing fraction of all avalanches 
and take the same values conjectured in Ref. 7. Be- 
yond these results, our analysis exemplifies a novel path 
towards the correct characterization of criticality which 
should be followed also for other models and phenomena 
out of equilibrium. 

Attempts to identify possible multiscaling features in 
avalanche distributions started soon after the introduc- 



Zj — Aji , where A is the discrete Laplacian. Thus, when 
toppling, site i loses four grains, giving one of them to 
each of its neighbors. At the border, dissipative Dirichlct 
conditions are assumed, so that one grain per toppling 
leaves the sandpile. Given a stable configuration (zk < Zc 
V k), the dynamics starts by adding a grain at a randomly 
chosen site Z (z/ — > z; -I- 1). U zi-\-l > Zc the site topples, 
and, if the case, other sites topple, as a consequence of 
this first instability, until a new stable configuration is 
reached. The sequence of all such topplings constitutes 
an avalanche. 

Grain addition is a slow driving mechanism compen- 
sated by border dissipation, which allows to reach the 
stationary SOC state after sufficiently many additions. 
If one considers a large number of bursts at stationar- 
ity, on average, one grain per event must be dissipated, 
because each avalanche originates from addition of one 
grain. However, this is still poor information on the 
dissipation mechanism. Indeed, outflow of grains oc- 
curs intermittently, concentrating on a small fraction of 
avalanches, which are separated by random sequences of 
non-dissipative bursts. This fraction approaches zero for 
L ^ OO. Thus, dissipating avalanches increase indefi- 
nitely with L their outflow. This fact has crucial conse- 
quences on the role played by this subset of avalanches 
on the global statistics. 

For the distribution P of the total number of topplings 
in an avalanche, s, FSS would imply an asymptotic form 
P{s,L) ~ s~^p[s/ L^), for s,L ^ OO, where D is a ca- 
pacity fractal dimension of the topplings. In order to 
better identify the role of rare, dissipating avalanches, we 
sampled only avalanches generated by grains added right 
at the center of the L x L box. In this way avalanches 
which dissipate a number of grains ti > at the borders 
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are always necessarily large, since they span a distance 
of order L. Such central seed setup is also most suit- 
able for our study of correlation functions. However as 
discussed below, our results can be obtained also with 
random additions over the whole box. Between different 
central seed additions we decorrelated the system with a 
sufficiently large number of arbitrary avalanches. 

We determined the outflow probability distribution 
0{n,L) for L = 32,64,128,256 and its moments 
J2nO{n,L)ni = ^Ez"f =< '^'^ >L, where N is the 
number of sampled avalanches (< 10*) and rii is the out- 
flow of the i-th one. Also the inverse first return fre- 
quency, i.e. the average number of avalanches, T, be- 
tween two successive nonzero outflows was determined. 
The intermittent character of the stationary state is best 
revealed by the law T{L) oc L^, which is satisfied with 
C — 0.50 ± 0.01. Intermittency and stationarity deter- 
mine a peculiar scaling of the moments of O. Indeed, the 
fraction of avalanches with n > is of order L^^. Thus, 
limq^o < >loc L^^. On the other hand, grain con- 
servation at stationarity also requires < n >]^= 1 — L'^ 
(one grain per avalanche added). Thus, by defining in 
general < >loc L^'^"'-i\ one expects tT„ > (o-„ < 0) 
for 0<g<l(q>l). We verified numerically that 
the average outflow, restricted to avalanches with n > 0, 
< n >iL= X^n >o number of avalanches 

with n > 0), grows as T, i.e. oc L''. These avalanches 
rarefy and grow in size for L oo. The most simple scal- 
ing to expect is < n"^ >l— L^^ < n'' >iloc L^^'^^^ for 
g > 0. This implies a constant gap cr„((7) — cr„ ((7— 1) = — C 
for q > 1, i.e. a linear behavior for cr„, consistent with 
FSS. Such a simple picture is confirmed by our numerical 
results. The behavior of cr„ is rather close to linear for 
all (7 > with a gap 2± —1/2 and ct„(0) ~ 1/2. Some 
previous determinations of in the literature fully agree 
with the present one |11 1^ and we conjecture ( = 1/2 
exactly. 

In order to elucidate the role of outflowing avalanches 
in determining the distribution P of the number of top- 
plings, we define gi{r,L) as the number of topplings in- 
duced at site r during the i-th avalanche. In view of the 
toppling rules, we clearly have: 



(1) 



where the Laplacian acts on r and the sum extends to the 
whole box. On the other hand, J2r 9ii^^ ^) ~ 
number of topplings in avalanche i. Thus, dimensional 
analysis alone would suggest: 



(2) 



i.e. cFsil) = ^2(7 + an{q). However, in spite of the strict 
linearity of cr„ assumed above, a constant gap behavior 
for as, as suggested by Eq. ||, is not acceptable. First of 
all, since s > O is not selective of outflowing avalanches, 
we can not have lim^^o < >lcx: L^^. Thus (Ts{0^) 7^ 



cr„(0+). Indeed, we know that avalanches with s > are 
a nonzero fraction of the total sampled for L ^ 00 
This implies limg_,o Cs (9) — 0. In addition, < s >l= 
J2r ^ di^' ^) >loc has been rigorously shown on the 
basis of stationarity and conservation 1). Thus, while 
Eq. § suggests as{q + 1) — (^siq) = —2.5 for all (7 > 0, 
basic properties of the non-equilibrium stationary state 
impose (Ts(1) — fTs(O) = —2 exactly. It remains to be 
decided whether a gap —2.5 still applies to at least part 
of the moments of P. Below we produce strong evidence 
in support of such a conclusion. 

Based on L = 64, 128, 256, 512, we extrapolated as as 
reported in Fig. |l|. We flnd limg_,o <^s{q) — 0, and, within 
an accuracy of lO"'^, crs(l) = —2, as expected. Most 
remarkably, moments of order <7 > 1 appear to conform 
very accurately to the constant gap —2.5 ± 0.1 expected 
on the basis of Eq. || (e.g., as{2) = 4.5 ±0.05). Moreover, 
for < (7 < 1, Us behaves nonlinear ly, clearly conflrming 
multiscaling for P. Outflowing avalanches alone obey 
simple scaling to very high accuracy. In Fig. ^ we report 
the moment exponent ais{q) for the distribution Pi of 
n > avalanches. Their dominance in the cumulative 
statistics of all avalanches is such that they impose their 
constant gap —2.5 as soon as (7 > 1. 
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FIG. 1. Plots of (Js (solid line) and ctis (dotted line); 
the dashed line has slope —2.50 ± 0.05 and y-intercept 
C, — 0.50 ± 0.01. Inset; plot of f{a) vs. a extrapolated from 
avalanche size distributions. Assuming f{a{oo)) = —1/2 we 
estimate a{oo) — —2.5 ± 0.1. 

An alternative way to discuss multiscaling is through 
the spectrum of singularity strengths /(a), i.e. the Leg- 
endre transform of as [p^ : a = das/dq, /(a) = —as + qa. 
a{q) and / are deflncd with reference to a saddle point 
(s = s*{q)) evaluation of < s'^ >l in the L ~> 00 limit 
{s*{q) ~ L-", P{s*,L) ~ L'^+f('^y). The inset in Fig. |l| 
shows the behavior of / as derived by directly transform- 
ing our L ~ 00 extrapolated data for as- Alternatively 
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one can extrapolate ensemble evaluations of / and a at 
finite L. The discrepancy of such different determina- 
tions allowed us to estimate their accuracy. The shape 
of /(a) is of course consistent with the above results for 
as'- in particular, the fact that /(— 5/2) — —1/2 is very 
well satisfied, agrees with the constant gap —2.5 expected 
for (Ts and with the fact that the linear continuation of 
the asymptotic as curve (dashed in Fig. |l|) intercepts the 
vertical axis a,t y — = 1/2. 

We verified that the above properties of P and as re- 
main substantially unaltered if sampled avalanches are 
created by random grain additions, occurring with uni- 
form probability on the whole box. In particular, n > 
is still a condition able to collect the subset of large 
avalanches dominating for g > 1 in the L — > oo limit. 
Of course, in this case the subset includes also many 
bursts of activity (e.g., boundary avalanches), whose im- 
portance becomes negligible in the thermodynamic limit. 
The central seed scheme offers the relevant advantage of 
more rapid convergence with comparable sizes and statis- 
tics. 

The constant gap ior q > 1 and the behavior of f{a) 
can be fully understood on the basis of the dominance 
of rare outflowing avalanches. In our sampling, such 
avalanches are also the largest, as far as distance spanned 
is concerned. We can write 

< iJ2 9ir, L)Y>L^^< {J2 9ir, L)f >ol (3) 

r r 

r 

where Nq is the number of not dissipating avalanches, 
and the subscripts of averages indicate restriction to the 
corresponding sets. In Eq. S of course, ^ ~ and 
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FIG. 2. Collapse fit for < g{r,L) >il. The inset re- 
ports a similar collapse for < g{r,L) >l- Data refer to 
L = 64, 128,256,512. 

For q ^ I, < g{v,L) >l and < g{r,L) >il in Eq. | 
satisfy very well ESS, but in different forms. Our data 
for < g{r, L) >l are consistent with a nearly logarithmic 
dependence on r/L (up to effects due to the upper and 
lower cutoffs). This is illustrated by the collapse plot in 
the inset of Fig. ||. Such logarithmic behavior oi < g >l 
agrees with the fact that this function has to coincide 
with the inverse Laplacian on the box. Indeed, due to 
the local conservation of grains, which holds on average 
for the whole sample of avalanches < <? >l must sat- 
isfy a Poisson equation with a unit source in the origin 
r = 0. On the other hand, for < g(r, L) >il one can not 
invoke local conservation. In fact, as shown by the col- 
lapse plot in Fig. 1^, our data are very well consistent with 
a form < g{r,L) >il— L^^'^gi{r/L), and the function gi 
is clearly not logarithmic. The peculiar anomalous scal- 
ing dimension 1/2 appearing in < g >il should be iden- 
tified with (. The collapses in Fig. ^ show that ti > 
avalanches are indeed much richer in topplings than the 
global average. The average number of waves in each out- 
flowing avalanche, coinciding with < g{0, L) >il is 
proportional to L^^^. In fact, if referred to all avalanches, 
such an average number grows only logarithmically with 
L 1^, as we verifled to high accuracy. Both terms on 
the r.h.s. of Eq. | give a contribution proportional to L 
when L —^ oo, as expected for the l.h.s.. 

The singularity strength /(a) is also fully consistent 
with the properties discussed above. /(— 2.5) = —1/2 
indicates that outflowing avalanches, which are a fraction 
7^-1/2 of the total, have a fractal dimension D = 2.5. 
In fact, within the multiscaling framework one obtains 
a whole continuum of fractal dimensions Dq — [(Ts(1) — 
as{q)]/[q— 1] for avalanches. Dq ranges in the whole 
interval (2, 2.5), for < g < 1, with Dq ^ 2 and Di = 2.5 
(Dq = for q > 1). 

Clearly, the very notion of standard scaling exponents 
can not adapt to P, due to the multiscaling concentrated 
at low q. However, in numerical work, P is usually ana- 
lyzed by assuming a ESS form, e.g. on the basis of col- 
lapse fits. Thus, it is important to ask what should be the 
ESS form of distribution which most accurately repro- 
duces the scaling of the moments of P. The optimal job 
in this respect is done by a distribution of the ESS form 
s~'^p{s/ L^), whose q-th moment exponent coincides with 
as{q) for q > 1, and with the function represented by the 
straight dashed line in Fig. ^ for q < 1. Of course, we 
put the fractal dimension equal to D — Di = 2.5 in 
p, because its opposite has to coincide with the asymp- 
totic slope of as- For such a scaling function one gets 
s^-Xs/L^) ~ L^(2-r) rpj^^g^ D[2-t) = 2 leads to 
T = 6/5. Some accurate numerical determinations based 
on ESS § are very close to this r, which remarkably co- 
incides with that conjectured in Ref. 7. The geometry of 
the various lines in Fig. |l]is very eloquent. In order to re- 
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duce the multiscaling of P to simple FSS, we have to shift 
to s = the contribution of n = avalanches in the his- 
togram of P. This amounts to bend the actual <Ts curve 
in the interval (0,1) into a straight segment (dashed), 
with intercept C on the vertical axis. Our results for Dq 
and T give an indication of the problems arising when, 
e.g., one tries to enforce FSS collapses of data for the 
simultaneous determination of D and r. Privileging col- 
lapse in the region of low s / produces a simultaneous 
lowering of both D and r, like when searching a solu- 
tion of Dq{2 — r) = 2 for low values of q. The extreme 
case is T = 1 with D = Dq = 2. Emphasizing collapse 
at high s /L^ produces higher determinations of both D 
and r. In this case a relatively poor sampling can lead to 
an overestimation of both exponents with respect to the 
above effective values. Indeed, we verified that and 
thus Dq tend to be overestimated systematically at large 

q- 

The standard exponents of rt > avalanches alone (en- 
semble 1) are n and D such that we can put Pi{s, L) — 
s~'^^Pi{s/L^). From our results and from the constant 
gap —2.5 for ais, it follows immediately D = 2.5 and 

n = 1. 

Summarizing, we showed that standard FSS does not 
hold for the ASM in 2D. Conservation of grains is 
guaranteed by an intermittent mechanism, with rare, 
large avalanches producing the outflowing current. These 
avalanches are such to fully determine the moments of P 
from the first one upwards. The multiscaling spectrum 
of singularity strengths for P has peculiar features as- 
sociated with the dominance of rare events. In particu- 
lar /(— 5/2) = —1/2 shows that rare avalanches have a 
fractal dimension equal to 5/2 and occur with frequency 
c>c L~^/^. These features are such to give r = 6/5 and 
D — 2.5 as effective exponents representing the multiscal- 
ing within an imposed simple FSS framework. Not sur- 
prisingly, FSS based numerical determinations of r and 
D are often quite close to the values mentioned above . 

In spite of the fact that the approach of Ref. 7 does 
not take into account multiscaling features, the values 
of r and D we propose coincide with those conjectured 
there. Not surprisingly, in view of our results, that con- 
jecture has met problems of numerical verification 
A revisitation of wave, or, rather, cluster [16| proper- 
ties in the light of the dominance of large avalanches and 
multiscaling can clarify the situation [ p^ . 

Crucial to the identification of our effective exponents 
are the constant gap — (2 + C) for q > 1 and the exact 
result ^^(l) = —2 Thus, it is precisely ^ = 1/2 that 
determines the conjectured t and D in the present frame- 
work. The critical state of the ASM in 2D is expected to 
correspond to that of a non-unitary conformal field the- 
ory with central charge c ~ —2 The compatibility of 
multiscaling with other non-unitary conformal field the- 
ories has been pointed out recently ||l^ . It is also known 
pgf that in a theory with c = —2 a correlator of disorder 
operators possesses exactly the scaling dimension 1/2. 

We expect intermittency and dominance of rare, large 



events to play a key role also in the physics of other SOC 
models. Indeed, ID sandpiles most often display pro- 
nounced multiscaling features and intermittency [pl]-p^ . 
Of course, the Abelian symmetry and the Laplacian char- 
acterization of the toppling dynamics made the anal- 
ysis of these features relatively easy here, to the ex- 
tent that exact properties of the model could be in- 
ferred. The case of ASM in higher dimensions, for which 
V < g{r, L) >loc seems to hold irrespective of d 
is also quite intriguing. 
Partial support from the European Network Contract 
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thank ISAS, Trieste, where part of this work was done. ^ 
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